
#include "heightfield.h"

#include <GL/glut.h>
#include "generador.h"
#include <cmath>

#include <cassert>

#include <limits>

// Se asume que existe una variable "numX" en el scope
#define ind(x, y) ((x) + ((y) * numX))

// Se asume que existe las variables "numX" y "z" en el scope
#define at(x, y) (z[((x) + ((y) * numX))])


void HeightField::add(int x, int y, float mx, float my) const {
   int indice = ind(x, y) * 3;

   float norma = sqrt(mx * mx + my * my + 1);
   normales[indice] = -mx/norma;
   normales[indice + 1] = -my/norma;
   normales[indice + 2] = 1.0/norma;
}

void HeightField::construirNormales() {
   /*
    * Se espera que el mapaZ este almacendado en columnas (for y: for x:)
    El mapa de normales sera guardado en el mismo orden.

    Sobre un (xy)[i], la normal sobre ese punto puede ser calculada.
    Sea la distribucion xy uniforme, esto es: 
    x[i] - x[i-1] = y[i] - y[i-1] = 1

    La pendiente mx de la recta sobre el eje x es (aproximadamente de orden 2):
    (z[i+1] - z[i-1]) / 2

    Idem para el eje y
    (z[j+1] - z[j-1]) / 2

    Donde z[i+1] es la componente z del punto x[i+1] y[i]
    Donde z[j+1] es en el punto x[i] y[i+1]

    Las rectas tangentes resultan:
    (x[i] y[i], z[i]) + (1, 0, mx) * u  (para el eje x)
    (x[i] y[i], z[i]) + (0, 1, my) * v  (para el eje y)

    El producto vectorial resulta ser un vector normal a la superficie
    i     j     x
    1     0     mx    ==> (-mx, -my, 1) Se retornanra una normal normalizada
    0     1     my

    En caso de no existir z[i+1] (un borde) usar
    z[i] - z[i-1] / 2       como pendiente en mx

    En otros casos:
    z[i+1] - z[i] / 2       como mx (falta z[i-1]
    z[j] - z[j-1] / 2       como my  (falta z[j+1]
    z[j+1] - z[j] / 2       como my  (falta z[j-1]
    */

   if(normales) delete[] normales;
   normales = new float[numX * numY * 3];

   float mx, my;

   //Esquina 0,0
   mx = (at(1, 0) - at(0, 0)) / 2.0;
   my = (at(0, 1) - at(0, 0)) / 2.0;

   add(0, 0, mx, my);

   //Borde 0,y
   for(int y = 1; y < numY - 1; ++y) {
      mx = (at(1, y) - at(0, y)) / 2.0;
      my = (at(0, y+1) - at(0, y-1)) / 2.0;

      add(0, y, mx, my);
   }

   //Esquina 0,Y
   int Y = numY - 1;
   mx = (at(1, Y) - at(0, Y)) / 2.0;
   my = (at(0, Y) - at(0, Y-1)) / 2.0;

   add(0, Y, mx, my);

   //El centro
   for(int x = 1; x < numX-1; ++x){

      //Borde x,0
      mx = (at(x+1, 0) - at(x-1, 0)) / 2.0;
      my = (at(x, 1) - at(x, 0)) / 2.0;

      add(x, 0, mx, my);

      for(int y = 1; y < numY-1; ++y) {

         //El centro, sin excepciones
         mx = (at(x+1, y) - at(x-1, y)) / 2.0;
         my = (at(x, y+1) - at(x, y-1)) / 2.0;

         add(x, y, mx, my);
      }

      //Borde x,Y
      mx = (at(x+1, Y) - at(x-1, Y)) / 2.0;
      my = (at(x, Y) - at(x, Y-1)) / 2.0;

      add(x, Y, mx, my);
   }

   //Esquina X,0
   int X = numX - 1;
   mx = (at(X, 0) - at(X-1, 0)) / 2.0;
   my = (at(X, 1) - at(X, 0)) / 2.0;

   add(X, 0, mx, my);

   //Borde X,y
   for(int y = 1; y < numY-1; ++y) {
      mx = (at(X, y) - at(X-1, y)) / 2.0;
      my = (at(X, y+1) - at(X, y-1)) / 2.0;

      add(X, y, mx, my);
   }

   //Esquina X,Y
   mx = (at(X, Y) - at(X-1, Y)) / 2.0;
   my = (at(X, Y) - at(X, Y-1)) / 2.0;

   add(X, Y, mx, my);

}

HeightField::HeightField(int numX, int numY, float densidad, 
      Generador &gen) : 
   numX(numX), 
   numY(numY),
   maxZ(std::numeric_limits<float>::min()),
   minZ(std::numeric_limits<float>::max()),
   z(0),
   densidad(densidad),
   normales(0),
   gen(gen) {
      construirZ();
   }

int HeightField::getX() const { return numX; }
int HeightField::getY() const { return numY; }

float HeightField::maxHeight() const { return maxZ; }
float HeightField::minHeight() const { return minZ; }

Coord HeightField::operator()(float u, float v) const {
   Coord coord(u*(numX-1), v*(numY-1), 0);

   int x = int(coord.x() + 0.5);
   int y = int(coord.y() + 0.5);

   if(x >= numX) { x = numX - 1; }
   if(y >= numY) { y = numY - 1; }

   coord.z() = at(x, y);

   return coord;
}

Coord HeightField::normal(float u, float v) const {
   int x = int(u*numX+ 0.5);
   int y = int(v*numY + 0.5);

   if(x >= numX) { x = numX - 1; }
   if(y >= numY) { y = numY - 1; }

   return Coord(&normales[ind(x, y) * 3]);
}

const float *HeightField::data() const {
   return z;
}

void HeightField::construirZ() {
   if(z) delete[] z;
   z = new float[numX * numY];

   int i = 0;
   for(int y = 0; y < numY; ++y) {
      float dy = y / densidad;
      for(int x = 0; x < numX; ++x) {
         assert(i < numX * numY);
         z[i] = gen(x / densidad, dy);
         maxZ = maxZ < z[i] ? z[i] : maxZ;
         minZ = minZ > z[i] ? z[i] : minZ;
         ++i;
      }
   }
   assert(i == numX * numY);

   construirNormales();
}


HeightField::~HeightField() {
   if(z) delete[] z;
   if(normales) delete[] normales;
}

